% center Aldh1a3 HCR signal in tdTomato cells from OTRcre x ai14
% resize to half resolution then save as .mat file with
%delta x y from original cell positions.

%% first run HCR_Cell_Localization_Batch_OTRcre_Slc6a5_Aldh1a3.m to create the HCR structure
%% Create data matrix, make box larger (center-80:center+80, then plot each one and select real center, record position).

path1='C:\Users\Tom�s\Documents\1. Regehr Lab\Data\HCR\OTRcrex Ai14 Slc6a5 Aldh1a3 12.7.20\Processing\4. position measurements\';
dir_to_search=[path1];
% pathtemp='C:\Users\Tom�s\Documents\1. Regehr Lab\Data\HCR\Registration and analysis\3.Channel colocalization\s1b2\';
path2='C:\Users\Tom�s\Documents\1. Regehr Lab\Data\HCR\OTRcrex Ai14 Slc6a5 Aldh1a3 12.7.20\Processing\1. Processed\';
path3='C:\Users\Tom�s\Documents\1. Regehr Lab\Data\HCR\'

% get the folder contents
d = dir(path1);
% remove all files (isdir property is 0)
dfolders = d([d(:).isdir]==1);
% remove '.' and '..'
dfolders = dfolders(~ismember({dfolders(:).name},{'.','..'}));

clear Cells ch1values ch2values ch3values

for j=1:numel(dfolders)

    Cells=HCR(j).CC(:,[1:2]);


ch1=imread([path2 dfolders(j).name '\' 'MAX_tdTomato_.tif'],1);% Fiji saved
ch2=imread([path2 dfolders(j).name '\' 'Slc6a5tdTomSubMatlabout.tiff'],1);% Matlab out
ch3=imread([path2 dfolders(j).name '\'  'Aldh1a3Matlabout.tiff'],1); %Matlab out

% figure
% imagesc(ch1)
% caxis([0 300])
% hold on
% scatter(CCs(:,1),CCs(:,2),10,'filled','r')

% ch1Mat=nan(121,121,size(Cells,1));
% ch2Mat=nan(121,121,size(Cells,1));
% ch3Mat=nan(121,121,size(Cells,1));
for i=2:size(Cells,1)
    ch1Mat(:,:,i)=ch1([round(Cells(i,2))-100:round(Cells(i,2))+100],round([round(Cells(i,1))-100:round(Cells(i,1))+100]));
    ch2Mat(:,:,i)=ch2([round(Cells(i,2))-100:round(Cells(i,2))+100],round([round(Cells(i,1))-100:round(Cells(i,1))+100]));
    ch3Mat(:,:,i)=ch3([round(Cells(i,2))-100:round(Cells(i,2))+100],round([round(Cells(i,1))-100:round(Cells(i,1))+100]));

end

if j==1
    cumch1Mat=ch1Mat;
    cumch2Mat=ch2Mat;
    cumch3Mat=ch3Mat;
else
end

% concatenate across slices
if ~(j==1)
    cumch1Mat=cat(3,cumch1Mat,ch1Mat);
    cumch2Mat=cat(3,cumch2Mat,ch2Mat);
    cumch3Mat=cat(3,cumch3Mat,ch3Mat);
else
end


end

% %% plot each cell HCR signal overlay using imshowpair
% figure
% for i=1001:size(cumch1Mat,3)
%     
%     subaxis(2,1,1)
%     
%   imshowpair(cumch1Mat(:,:,i),cumch3Mat(:,:,i));
%   vline(81,'w')
%   hline(81,'w')
%    title(num2str(i))
%     subaxis(2,1,2)
%     imagesc(cumch3Mat(:,:,i))
%      colormap([cmap('magenta',100,0,50)]);
%    caxis([200 1500])
%      vline(81,'w')
%   hline(81,'w')
%      axis square
%     [xi(i),yi(i)] = getpts
% 
%     pause
% end
% %save progress lol
% %I have done until 517
% saveX=xi;
% saveY=yi;
% 
% deltaX=80-xi; % xi yi were taken from 80 x 80 matrix. delta should be useful for any matrix.
% deltaY=80-yi;

% x100=xi+20;
% y100=xi+20;

if sum(size(Cells)==size(HCR(j).CC(:,[1:2])))==2  % just means if i'm analyzing CCs that were recentered, load the delta x and delta y
  load('C:\Users\Tom�s\Documents\1. Regehr Lab\Data\HCR\OTRcrex Ai14 Slc6a5 Aldh1a3 12.7.20\Processing\4. position measurements\deltaXdeltaYAldh1a3.mat')
  deltaX(deltaX>99)=0; % don't shift cells that have wrong delta X or delta Y (only 2 cells, I think). This is a quick fix but not ideal.
  deltaY(deltaY>99)=0;
 
  x100=deltaX+100;
y100=deltaY+100;
end

%%Trim
for i=1:size(x100,2)
    ch1MatTrim(:,:,i)=cumch1Mat([round(x100(i))-60:round(x100(i))+60],[abs(round(y100(i)))-60:abs(round(y100(i)))+60],i);
    ch1MatTR(:,:,i)=imresize(ch1MatTrim(:,:,i),[61 61],'bicubic');
    ch2MatTrim(:,:,i)=cumch2Mat([round(x100(i))-60:round(x100(i))+60],[abs(round(y100(i)))-60:abs(round(y100(i)))+60],i);
     ch2MatTR(:,:,i)=imresize(ch2MatTrim(:,:,i),[61 61],'bicubic');
    ch3MatTrim(:,:,i)=cumch3Mat([round(x100(i))-60:round(x100(i))+60],[abs(round(y100(i)))-60:abs(round(y100(i)))+60],i);
     ch3MatTR(:,:,i)=imresize(ch3MatTrim(:,:,i),[61 61],'bicubic');

end



% save([path1 'deltaxdeltayAldh1a3.mat'],'deltaX', 'deltaY')

% %% new section for removing slc artifact cells.
% figure
% for i=1:size(ch2MatTR,3)
%     imagesc(ch2MatTrim(:,:,i))
%     title(num2str(i))
%     caxis([200 1500])
%     colormap([cmap('cyan',100,0,0)]);
%    
%     prompt=' 1 good, 0 bad' 
%  bad(i)=input(prompt)
% end
%    
   
   %% plot
   
   radius=32;
ylimitprofiles=600;
   
   figure
c1Ax=subaxis(2,3,1)
  imagesc(nanmean(ch1MatTR,3))
  caxis([0 12000])
    y1=colormap(c1Ax,[cmap('gray',100,0,20)]);

  axis square
  axis off
  subaxis(2,3,4)
%   plot(mean(mean(cumch1Mat,3),1))
meanGreen=nanmean(ch1MatTR,3);
greenHor=mean(meanGreen([29:33],:),1);
greenVert=mean(meanGreen(:,[29:33]),2);
greenProf=mean([greenVert';greenHor],1)
  hold on
  plot(greenProf,'k')
%   plot(mean(mean(cumch1Mat,3),2)) %average profiles in x and y to decide on boundary of circular mask
  
  ylim([0 ylimitprofiles*40/2])
%   vline([30-radius 30+radius],'--k')
  %draw axis for visual purpose
  vline(0,'k')
  hline(0,'k')
  axis square
  axis off
  
 c2Ax= subaxis(2,3,2)
  imagesc(nanmean(ch2MatTR,3))
  caxis([0 350])
  y2=colormap(c2Ax,[cmap('cyan',100,0,50)]);
  axis square
  axis off
  subaxis(2,3,5)
%   plot(mean(mean(cumch2Mat,3),1))
%   hold on
%   plot(mean(mean(cumch2Mat,3),2)) %average profiles in x and y to decide on boundary of circular mask
meanCyan=nanmean(ch2MatTR,3);
CyanHor=mean(meanCyan([29:33],:),1);
CyanVert=mean(meanCyan(:,[29:33]),2);
CyanProf=mean([CyanVert';CyanHor],1)
  hold on
  plot(CyanProf,'k')
   ylim([0 ylimitprofiles/2])
%    vline([30-radius 30+radius],'--k')
   %draw axis for visual purpose
  vline(0,'k')
  hline(0,'k')
  axis square
  axis off
  
   c3Ax= subaxis(2,3,3)
  imagesc(nanmean(ch3MatTR,3))
  
  caxis([0 350])
    y3=colormap(c3Ax,[cmap('magenta',100,0,50)]);

  axis square
  axis off
  subaxis(2,3,6)
%   plot(mean(mean(cumch3Mat,3),1))
%   hold on
%   plot(mean(mean(cumch3Mat,3),2)) %average profiles in x and y to decide on boundary of circular mask
meanMagenta=nanmean(ch3MatTR,3);
MagentaHor=mean(meanMagenta([29:33],:),1);
MagentaVert=mean(meanMagenta(:,[29:33]),2);
MagentaProf=mean([MagentaVert';MagentaHor],1)
  hold on
  plot(MagentaProf,'k')
   ylim([0 ylimitprofiles/2])
%    vline([30-radius 30+radius],'--k')
   %draw axis for visual purpose
  vline(0,'k')
  hline(0,'k')
  axis square
  axis off

  
  
  %%   %% masking
    
% Create a logical image of a circle with specified
% diameter, center, and image size.
% First create the image.
imageSizeX = 61;
imageSizeY = 61;
[columnsInImage, rowsInImage] = meshgrid(1:imageSizeX, 1:imageSizeY);
% Next create the circle in the image.
centerX = ceil(imageSizeX / 2); % in the middle.
centerY = ceil(imageSizeY / 2);
radius = 50; %based on profiles

circlePixelsLogical = (rowsInImage - centerY).^2 ...
    + (columnsInImage - centerX).^2 <= radius.^2;
% circlePixels is a 2D "logical" array.
% Now, display it.
figure
image(circlePixelsLogical) ;
colormap([0 0 0; 1 1 1]);
title('Binary image of a circle');
axis square;

circlePixels01=double(circlePixelsLogical);
test=circlePixels01.*double(ch1MatTR); %mask out crap outisde cells 
figure
imagesc(nanmean(test,3))
caxis([0 50])
axis square



%% take averages of each cell for each channel inside the circular mask
clear ch1values
clear ch2values
clear ch3values
clear ch1tsne
clear ch2tsne
clear ch3tsne
for c=1:size(ch1MatTR,3)
    tempdata1=ch1MatTR(:,:,c);
%     ch1tsne(c,:)=tempdata1(circlePixelsLogical);
    ch1values(c)=mean(tempdata1(circlePixelsLogical));
    
     tempdata2=ch2MatTR(:,:,c);
%       ch2tsne(c,:)=tempdata2(circlePixelsLogical);
    ch2values(c)=mean(tempdata2(circlePixelsLogical));
    
     tempdata3=ch3MatTR(:,:,c);
%       ch3tsne(c,:)=tempdata3(circlePixelsLogical);
    ch3values(c)=mean(tempdata3(circlePixelsLogical));
    
end

figure
limVal=500;
subaxis(1,3,1)
scatter(ch1values,ch2values,10,'filled','k')
xlabel('tdTomato')
ylabel('Slc6a5')
axis square
xlim([0 30e3])
ylim([0 limVal])
subaxis(1,3,2)
scatter(ch1values,ch3values,10,'filled','k')
xlabel('tdTomato')
ylabel('Aldh1a3')
xlim([0 30e3])
ylim([0 limVal])
axis square
subaxis(1,3,3)
scatter(ch2values,ch3values,10,'filled','k')
xlabel('Slc6a5')
ylabel('Aldh1a3')
xlim([0 limVal])
ylim([0 limVal])
axis square

%% new section to try out the combination of circular mask with threshold mask based on tdTomato fluorescence


for i=1:size(ch1MatTR,3)
a=ch1MatTR(:,:,i);
back=prctile(a,5);
maxCelle=max(max(a([22:37],[22:37])));
b=a>back*3 & a<maxCelle;
threshch1Masks(:,:,i)=b & circlePixelsLogical;

% imshowpair(double(b),a,'montage')
% pause
end

% take mean values
for c=1:size(ch1MatTR,3)
    tempdata1=ch1MatTR(:,:,c);
%     ch1tsne(c,:)=tempdata1(circlePixelsLogical);
    ch1valuesT(c)=mean(tempdata1(threshch1Masks(:,:,c)));
    
     tempdata2=ch2MatTR(:,:,c);
%       ch2tsne(c,:)=tempdata2(circlePixelsLogical);
    ch2valuesT(c)=mean(tempdata2(threshch1Masks(:,:,c)));
    
     tempdata3=ch3MatTR(:,:,c);
%       ch3tsne(c,:)=tempdata3(circlePixelsLogical);
    ch3valuesT(c)=mean(tempdata3(threshch1Masks(:,:,c)));
    masksize(c)=sum(sum(threshch1Masks(:,:,c)));
    
end
%% scatter plots with histograms below
figure
limVal=500;
tdTomlog=ch1values>0;
subaxis(2,2,1)
scatter(ch1valuesT(tdTomlog),ch2valuesT(tdTomlog),10,'filled','k','MarkerFaceAlpha',.7,'MarkerEdgeAlpha',.7)
xlabel('tdTomato')
ylabel('Slc6a5')
axis square
xlim([0 30e3])
ylim([0 limVal])
subaxis(2,2,2)
scatter(ch1valuesT(tdTomlog),ch3valuesT(tdTomlog),10,'filled','k','MarkerFaceAlpha',.7,'MarkerEdgeAlpha',.7)
xlabel('tdTomato')
ylabel('Aldh1a3')
xlim([0 30e3])
ylim([0 limVal])
axis square

subaxis(2,2,3) % histogram of number of cells vs tdTomato, binned.
histogram(ch1valuesT(tdTomlog))
xlim([0 30e3])
axis square
box off
% Fixed tdTom binning has issues:
% subaxis(2,2,4) % average slc & Aldh1a3 per tdTom bin
% binsize=1000;  
% x=binsize/2:binsize/2:3e4-binsize/2;
% binEdgesL=x-binsize/2;
% binEdgesR=x+binsize/2;
% meanSlc6a5tdTbin=[];
% meanAldh1a3tdTbin=[];
% for i=1:numel(x)
%     templog=ch1valuesT>binEdgesL(i) & ch1valuesT<binEdgesR(i);
%     i
%     sum(templog)
%     meanSlc6a5tdTbin(i)=nanmedian(ch2valuesT(templog));
%     meanAldh1a3tdTbin(i)=nanmedian(ch3valuesT(templog));
% end
% plot(x,medfilt1(meanSlc6a5tdTbin,1,'omitnan'),'linewidth',2,'color',[21 231 255]./255)
% hold on
% plot(x,medfilt1(meanAldh1a3tdTbin,1,'omitnan'),'linewidth',2,'color',[206 0 244]./255)
% axis square
% box off

% new version of average signal per bin but using bins of number of cells
% instead of ranges of tdTom values. This solves issues with so many cells
% in low values and few cells in higher tdTom. so it's a sort of rolling
% binning where X value is the mean tdTom in that bin instead of a fixed
% value.


subaxis(2,2,4) % average slc & Aldh1a3 per tdTom bin (100 cells)

binsizeCells=50; %number of cells per bin.
[b,idx]=sort(ch1valuesT); % sort tdTomato values.
nBins=floor(numel(ch1valuesT)/binsizeCells);

meanSlc6a5tdTbin=[];
meanAldh1a3tdTbin=[];
x=[];
for i=1:nBins
    
    if ~(i==nBins)
    tempIdx=idx([1:binsizeCells]+binsizeCells*(i-1));
    else
        tempIdx=idx(1+binsizeCells*(i-1):end);
    end
    x(i)=nanmean(ch1valuesT(tempIdx))
    meanSlc6a5tdTbin(i)=nanmedian(ch2valuesT(tempIdx));
    meanAldh1a3tdTbin(i)=nanmedian(ch3valuesT(tempIdx));
end
plot(x,medfilt1(meanSlc6a5tdTbin,1,'omitnan'),'linewidth',2,'color',[21 231 255]./255)
hold on
plot(x,medfilt1(meanAldh1a3tdTbin,1,'omitnan'),'linewidth',2,'color',[206 0 244]./255)
xlim([0 3e4])
axis square
box off
% title(['binsize= ' num2str(binsizeCells) ' Cells'])






%% Remake average HCR plot with cells above a certain tdTomato threshold

%%
% make nice figure 
%display all the cells stacked together to look at data distribution

tdTomlog=ch1values>0;
ylimitprofiles=300;
figure
c1Ax=subaxis(2,3,1)
  imagesc(nanmean(ch1MatTR(:,:,tdTomlog),3))
  caxis([0 20000])
    y1=colormap(c1Ax,[cmap('gray',100,0,0)]);

  axis square
  axis off
  subaxis(2,3,4)
%   plot(mean(mean(cumch1Mat,3),1))
meanGreen=nanmean(ch1MatTR(:,:,tdTomlog),3);
greenHor=nanmean(meanGreen([29:33],:),1);
greenVert=nanmean(meanGreen(:,[29:33]),2);
greenProf=nanmean([greenVert';greenHor],1)
  hold on
  plot(greenProf,'k')
%   plot(mean(mean(cumch1Mat,3),2)) %average profiles in x and y to decide on boundary of circular mask
  
%   ylim([0 ylimitprofiles*40/2])
%   vline([30-radius 30+radius],'--k')
  %draw axis for visual purpose
  ylim([0 20000])
  vline(0,'k')
  hline(0,'k')
  axis square
  axis off
  
 c2Ax= subaxis(2,3,2)
  imagesc(nanmean(ch2MatTR(:,:,tdTomlog),3))
  caxis([0 ylimitprofiles])
  y2=colormap(c2Ax,[cmap('cyan',100,0,50)]);
  axis square
  axis off
  subaxis(2,3,5)
%   plot(mean(mean(cumch2Mat,3),1))
%   hold on
%   plot(mean(mean(cumch2Mat,3),2)) %average profiles in x and y to decide on boundary of circular mask
meanCyan=nanmean(ch2MatTR(:,:,tdTomlog),3);
CyanHor=nanmean(meanCyan([29:33],:),1);
CyanVert=nanmean(meanCyan(:,[29:33]),2);
CyanProf=nanmean([CyanVert';CyanHor],1)
  hold on
  plot(CyanProf,'k')
   ylim([0 ylimitprofiles])
%    vline([30-radius 30+radius],'--k')
   %draw axis for visual purpose
  vline(0,'k')
  hline(0,'k')
  axis square
  axis off
  
   c3Ax= subaxis(2,3,3)
  imagesc(nanmean(ch3MatTR(:,:,tdTomlog),3))
  
  caxis([0 ylimitprofiles])
    y3=colormap(c3Ax,[cmap('magenta',100,0,50)]);

  axis square
  axis off
  subaxis(2,3,6)
%   plot(mean(mean(cumch3Mat,3),1))
%   hold on
%   plot(mean(mean(cumch3Mat,3),2)) %average profiles in x and y to decide on boundary of circular mask
meanMagenta=nanmean(ch3MatTR(:,:,tdTomlog),3);
MagentaHor=mean(meanMagenta([29:35],:),1);
MagentaVert=mean(meanMagenta(:,[28:34]),2);
MagentaProf=mean([MagentaVert';MagentaHor],1)
  hold on
  plot(MagentaProf,'k')
   ylim([0 ylimitprofiles])
%    vline([30-radius 30+radius],'--k')
   %draw axis for visual purpose
  vline(0,'k')
  hline(0,'k')
  axis square
  axis off



  %% Make cumulative plots for different tdTomato thresholds
  
  figure;histogram(ch2valuesT(ch1values>0),'DisplayStyle','stairs','Normalization','cdf','binwidth',1,'EdgeColor',[21 231 255]./255);
  hold on
  histogram(ch3valuesT(ch1values>0),'DisplayStyle','stairs','Normalization','cdf','binwidth',1,'EdgeColor',[206 0 244]./255);
  
  
